Beta Regression

Code
library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)

theme_set(theme_pubr(legend = "bottom"))


nnns_imputed<- readRDS("../Haojia_work/nnns_imputed.rds")

read imputed data

Code
dat=complete(nnns_imputed,1)
attach(dat);Percent_of_feeds_taken_by_mouth_at_discharge
  [1] 1.00 0.19 0.00 0.01 0.00 0.00 0.06 0.00 0.00 0.00 0.50 0.00 0.00 1.00 0.50
 [16] 0.02 0.00 0.00 0.29 0.00 0.00 0.01 0.14 0.30 0.06 0.41 0.95 0.01 0.17 1.00
 [31] 0.34 0.50 0.00 0.00 0.49 0.00 0.49 0.00 0.83 0.03 0.00 0.00 0.02 0.00 0.00
 [46] 0.00 1.00 1.00 0.00 0.00 0.20 0.01 0.00 1.00 0.05 0.00 0.14 0.00 0.12 1.00
 [61] 0.80 0.00 0.90 0.25 0.05 0.02 0.00 0.25 0.05 0.00 0.00 0.15 0.00 0.00 0.08
 [76] 0.00 0.03 0.30 0.23 0.01 0.50 0.09 0.10 1.00 0.44 0.80 0.50 0.22 0.15 0.00
 [91] 0.00 0.50 0.30 0.33 0.20 0.42 0.00 0.03 0.12 0.06 0.67 0.00 0.00 0.00 0.14
[106] 0.00 0.11 0.13 1.00 1.00 0.09 0.50 0.00 1.00 0.00

Use ZOIB package

ZOIB Model

\[ f(y_i|\eta_i) = \begin{cases} p_i, & y_i =0 \\ (1-p_i)q_i, & y_i =1 \\ (1 − p_i)(1 − q_i)Beta(\alpha_{i_1}, \alpha_{i_2}) &y_i \in (0, 1) \end{cases} \]

\[ \begin{aligned} logit\left(\mu^{(0,1)} = \frac{\alpha_1}{\alpha_1+ \alpha_2}\right) = \mathbf x_1 \boldsymbol \beta_1 + I_1(\mathbf z_{1,i}\gamma_{1,i}) \\ log\left(v_i = \alpha_1+ \alpha_2\right) = \mathbf x_2 \boldsymbol \beta_2 + I_2(\mathbf z_{2,i}\gamma_{2,i}) \\ logit\left(p_i\right) = \mathbf x_3 \boldsymbol \beta_3 + I_1(\mathbf z_{3,i}\gamma_{3,i}) \\ logit\left(q_i\right) = \mathbf x_4 \boldsymbol \beta_4 + I_1(\mathbf z_{4,i}\gamma_{4 ,i}) \end{aligned} \]

Fit ZOIB model

Code
set.seed(11282023) #Set seed, bayesian model!

tictoc::tic()
pre_op_model =zoib(Percent_of_feeds_taken_by_mouth_at_discharge #Response
            ~Pre_Op_NNNS_attention_score+
              Length_of_intubation_days+
              Genetic_Syndrome_or_Chromosomal_Abnormality+
              Cardiac_Anatomy+
              Age_at_Surgery_days+
              Female
              #x1 design matrix
       | 1 #x2 design matrix- estimating variance
     |Pre_Op_NNNS_attention_score+
              Length_of_intubation_days+
              Genetic_Syndrome_or_Chromosomal_Abnormality+
              Cardiac_Anatomy+
              Age_at_Surgery_days+
              Female #x3 design matrix
     |Pre_Op_NNNS_attention_score+
              Length_of_intubation_days+
              Genetic_Syndrome_or_Chromosomal_Abnormality+
              Cardiac_Anatomy+
              Age_at_Surgery_days+
              Female, #x4 design matrix
     data = dat,
     n.response = 1,
     zero.inflation = TRUE,
     one.inflation = TRUE,
     link.mu = "logit",
     link.x0 = "logit",
     link.x1 =  "logit",
     random = 0,
     n.chain = 4, # number of MCMC chains
     n.iter = 5000, #number of iterations before burn/thin
     n.thin = 2, # thinning period
     n.burn = 200, # burn-in period 
     seeds = c(11,29,20,23) #vector of seeds for chains to make reproducable
     )
tictoc::toc()
saveRDS(pre_op_model, file="model1.RData")

Interpreting output:

b: vector of estimates from Eqn 1; that is, g(mu) = xb*b +z*gamma

d: vector of estimates from Eqn 2; that is, log(eta) = xd*d+z*gamma

b0: vector of estimates from Eqn 3; that is, g(p0) = x0*b0 +z*gamma

b1: vector of estimates from Eqn ;4 that is, g(p1) = x1*b1+z*gamma

Code
model1 =readRDS(file="model1.RData")

model1$coeff |> summary()

Iterations = 1:2400
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 2400 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean       SD  Naive SE Time-series SE
b[1]   -0.40816  0.74818 0.0076361      0.0075400
b[2]    0.08670  0.16883 0.0017231      0.0017843
b[3]   -0.12256  0.06370 0.0006502      0.0007321
b[4]    0.01044  0.38861 0.0039662      0.0043701
b[5]    0.51766  0.39137 0.0039944      0.0040275
b[6]   -0.22897  0.35027 0.0035749      0.0040599
b[7]   -0.01677  0.03065 0.0003128      0.0003285
b[8]   -0.44165  0.30822 0.0031458      0.0033937
b0[1]  -1.88085  1.22359 0.0124882      0.0130445
b0[2]  -0.14518  0.25755 0.0026286      0.0026731
b0[3]   0.27908  0.09623 0.0009821      0.0010380
b0[4]   0.31831  0.59714 0.0060946      0.0063665
b0[5]   1.54558  0.57721 0.0058911      0.0065491
b0[6]   1.07857  0.58879 0.0060093      0.0069353
b0[7]  -0.04440  0.05236 0.0005344      0.0005509
b0[8]  -0.70889  0.49425 0.0050444      0.0052464
b1[1]  -5.26454  2.52180 0.0257380      0.0299346
b1[2]   0.60754  0.55635 0.0056782      0.0067350
b1[3]  -0.38975  0.21021 0.0021454      0.0031057
b1[4] -72.91600 41.62390 0.4248222      7.9233266
b1[5]  -0.09163  1.61224 0.0164549      0.0211416
b1[6]  -0.49582  1.00420 0.0102491      0.0123597
b1[7]   0.24888  0.10494 0.0010711      0.0013713
b1[8]   1.47877  0.88376 0.0090198      0.0100330
d       1.01102  0.17214 0.0017569      0.0019064

2. Quantiles for each variable:

            2.5%        25%       50%        75%      97.5%
b[1]    -1.85699   -0.90675  -0.40844   0.087590  1.0715096
b[2]    -0.23977   -0.02749   0.08513   0.200625  0.4193043
b[3]    -0.25062   -0.16489  -0.12183  -0.079148 -0.0006274
b[4]    -0.75772   -0.24505   0.01394   0.272360  0.7683315
b[5]    -0.25352    0.26030   0.52092   0.780517  1.2806785
b[6]    -0.91720   -0.46624  -0.22671   0.007560  0.4539170
b[7]    -0.07966   -0.03713  -0.01625   0.003648  0.0429441
b[8]    -1.05363   -0.64703  -0.44025  -0.230532  0.1575015
b0[1]   -4.34565   -2.68888  -1.88189  -1.035950  0.4109481
b0[2]   -0.65047   -0.31737  -0.14961   0.029522  0.3649150
b0[3]    0.10721    0.21125   0.27318   0.339820  0.4848661
b0[4]   -0.84659   -0.07684   0.31247   0.717412  1.4897533
b0[5]    0.43129    1.15243   1.53832   1.928722  2.6996217
b0[6]   -0.05431    0.68697   1.08013   1.465759  2.2500901
b0[7]   -0.15423   -0.07724  -0.04288  -0.008099  0.0526682
b0[8]   -1.68330   -1.04233  -0.70207  -0.367095  0.2396023
b1[1]  -10.69062   -6.84610  -5.12750  -3.502697 -0.7335187
b1[2]   -0.43126    0.23030   0.59253   0.959609  1.7767425
b1[3]   -0.84252   -0.52501  -0.37801  -0.241454 -0.0087959
b1[4] -150.35996 -106.74572 -73.38067 -37.844099 -7.5182184
b1[5]   -3.60955   -1.03421   0.03328   0.992953  2.7713079
b1[6]   -2.49611   -1.15544  -0.49286   0.190776  1.4089553
b1[7]    0.06095    0.17607   0.24163   0.315919  0.4728664
b1[8]   -0.22888    0.88836   1.45247   2.053894  3.2583483
d        0.66224    0.89529   1.01546   1.129489  1.3390093

Check convergence

Traceplots

Code
traceplot(model1$coeff)

Autocorrelation plots

Code
autocorr.plot(model1$coeff)